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The stable structures and melting properties of ion clouds in isotropic octupole traps are inves- 
tigated using a combination of semi-analytical and numerical models, with a particular emphasis 
at finite size scaling effects. Small-size clouds are found to be hollow and arranged in shells cor- 
responding approximately to the solutions of the Thomson problem. The shell structure is lost in 
clusters containing more than a few thousands of ions, the inner parts of the cloud becoming soft 
and amorphous. While melting is triggered in the core shells, the melting temperature unexpectedly 
follows the rule expected for three-dimensional dense particles, with a depression scaling linearly 
with the inverse radius. 
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I. INTRODUCTION 

The storage of cold ions in electrostatic and magnetic 
fields (Penning traps) or in radio-frequency electric fields 
(Paul traps) has becomepossible with the advent of laser 
cooling techniques [E H IS El- For few-particle sys- 
tems in quadrupole traps, this has lead to many appli- 
cations ranging from high resolution spectroscopy and 
optical frequency standards 0, Q to quantum informa- 
tion [3, [a Q and tests on the possible variations of the 
fundamental constants [ToL [ll| . While these effective 
harmonic traps allow focussing the ions at the center, 
a greater number of ions with reduced rf driven motion 
can be stored in higher-order confinements, with possible 
uses as microwave clocks for deep space navigation (l2j 
or to control cold chemical reactions [13[ of astrophysical 
relevance. 

Large samples of ions also offer a practical realiza- 
tion of the classical many-body Coulomb problem with 
interesting collective properties such as phase transi- 
tions. Crystallization of ion clouds has been observed in 
quadrupole traps [3, [TH, Ojj], in Penning traps [j, 17 1, 
and theoretically studied by several groups [I3 [l8 j. 



There, ion clouds crystallize into well-defined shells at 
small sizes 18], and then into the bcc Wigner crystal in 
samples exceeding about 10 4 ions Melting can be 

experimentally triggered by varying the trap parameters 
[20(, and proceeds by separate radial and orientational 
mechanisms in small clusters 1211. Schiffer [23 has re- 



ported from molecular dynamics (MD) simulations that 
the large clouds melt from the surface, resulting in a melt- 
ing temperature that varies linearly with inverse cluster 
radius. 

In comparison, the properties of cold ion clouds in 
higher-order traps are much less documented, despite 
specific investigations for octupole traps [23|, Hj] . Using 
various semi-analytical and numerical methods, we show 
in this article that the stable structures of ion clouds in 
octupole traps are generally made of multiple distinct 
shells only when the number of ions does not exceed a 
few thousands. Above this approximate size, layering is 
progressively lost except in the outermost regions of the 
cloud. We also found that these clusters melt from the 
core and, quite unexpectedly, exhibit a depression in the 
melting point that scales linearly with the cloud radius. 

In the next section, the model is described and the 
structural properties of clusters are investigated. In 
Sec. IIII1 the finite temperature aspects are covered, and 
some concluding remarks are finally given in Sec. IIV1 



II. MODEL AND STATIC PROPERTIES 

The system we are investigating consists of N identical 
ions with charge q and mass m, trapped in an isotropic 
octopole trap. We denote by f2 the radiofrequency of the 
electric field and by E (r) its amplitude at position r. We 
assume that the adiabatic approximation holds and that 
the macromotion is driven by the so-called pseudopoten- 
tial $t(t) as 
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We further assume that the rf-driven motion that su- 
perposes to the macromotion can be neglected for the 
present purposes. For simplicity of the following analy- 
sis, the electric field derives from a purely radial three- 
dimensional octupole potential £ — — gradV^ with V = 
Vo x r 4 , r being the distance from the trap center and Vq 
a constant. The pseudopotential then scales as r 6 and 
the potential energy felt by an assembly of ions can thus 
be written as 



*({*}) = *T({r i }) + *c({r i }), 



(2) 



i<j 



where we have denoted by and the distances of ion 
i to the center of the trap and to particle j, respectively. 
The constants A and B in the previous equation can be 
removed by scaling of the quantities r — * r = and <I> — > 
<3> = $/jB, with 7 = (A/B) 1 / 7 . In the following, reduced 
units are thus chosen for both distances and energies, 
which amounts to using A = B = 1 in Eq. @. 

We have first located the structures that globally min- 
imize Eq. ^ at fixed size N < 2 00, employing the basin- 
hopping Monte Carlo method [25[ successfully used in 
previous related work [26| . The ions arrange into a single 
spherical thin shell for N < 109, and two shells above this 
size. This shell structure can be fruitfully used to predict 
the minima for larger clouds, assuming that the radial 
density is written as a sum over M concentric shells of 
zero thickness. Such a shell model was initially developed 
for the quadrupole trap [2^ | and improved for intra-shell 
correlations [281 ]. The energy to be minimized, $ s hcii, 
is a function of the radii {Ri} of the shells, each shell i 
carrying Ni ions (plus one possible ion at the center, No): 
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In the above equation, we have denoted by E 1 ( ^ tra, (N) the 
intra-shell Coulomb energy in which the particles lie on 
the unit sphere. This energy is minimized at the solutions 
of the Thomson proble m [2911 , which have been tabulated 
up to rather large sizes [301 ] 



E^ a (N) = E Th (N). 



(4) 



Alternatively, for large numbers of ions, an asymptotic 
expression can be substituted for the intra-shell Coulomb 
energy [3l[ as 



(5) 



where the a\N contribution accounts for intra-shell 
correlations (neglecting correlations in a mean-field ap- 
proach would yield N(N — l)/2). Following the recent 
results of Cioslowski and Grzebielucha [3l[ the parameter 
a was taken as 1.10610 to account for correlations. The 



above energy <E> s h e ii({-Ri, Ni}) can be exactly minimized, 
resulting in the expression for the radius of shell i 



Ri — 



J5j? tra (JVi) 
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(6) 



The remaining minimization of $ s heii with respect to the 
{Ni} must then be carried out numerically under the 
constraint of a fixed total number of ions. 

When the tabulated optimal Thomson energies are 
used in place for E™ tr& (N), the minimization problem 
is variational and provides rigorous upper bounds to the 
exact values. If the asymptotic expression of Eq. ([5]) is 
employed instead, energies lower than the numerically 
exact minima may be reached due to the approximate 
nature of this asymptotic form. The optimized ener- 
gies, outer radius, and shell arrangements obtained from 
Monte Carlo minimization are listed in Table [I] for se- 
lected cluster sizes. In this table, the predictions of the 
shell models in which the intra-shell energies are either 
taken from the tabulated Thomson minima [30j (discrete 
model) or from the asymptotic expression (continuous 
model) are also given. 

Energies, shell radius and ion arrangements agree very 
well between the three methods, the shell models some- 
times producing differences of ±1 ion in the shells for the 
larger clusters containing 180 ions or more. Energies and 
radii obtained from the tabulated Thomson minima are 
in very good agreement with the exact results until two 
shells are formed above 100 ions. The agreement is not as 
good when the asymptotic form for the Thomson energies 
is used, the energy being slightly underestimated. This 
underestimation suggests that the value of the parameter 
a used in the shell model, which was taken from extrap- 
olations of the Thomson model to the large sizes regime 
[3l|, may be slightly excessive at small sizes. The greater 
validity of the current value of a to large sizes is also 
indicated by the relative error between the shell energy 
and the Monte Carlo data, which for the systems consid- 
ered in Table [J decreases from 0.5% to less than 0.01% as 
the number of ions increases from 10 to 200. Despite this 
systematic error, we generally find that the exact energy 
obtained by minimization lies inbetween the predictions 
of the two shell models, and that the outer radius is cor- 
rectly reproduced (within 0.2%) by both models. 

At this point, it is important to stress that the agree- 
ment between the shell models and global optimization 
is essentially due to the correct account of correlations, 
which are implicitly included in the Thomson energies 
ETh{N) or explicitly in the asymptotic expression of 
Eq ([5]) through the term SN = ay/~N. For comparison, 
the mean-field treatment predicts that the structure of 
the 100-ion cluster would have 9 shells. 

The asymptotic shell model can be further exploited 
in the larger sizes regime 10 3 < N < 10 5 , where the op- 
timal Thomson energies are not systematically available. 
The shell radii predicted by this model are represented 



3 





Minimization 


Thomson shell model 


Asymptotic shell model 


Size 


Energy 


Outer radius 


Arrangement 


Energy 


Outer radius 


Arrangement 


Energy 


Outer radius 


Arrangement 


10 


41.624 


0.917 


(10) 


41.624 


0.917 


(10) 


41.399 


0.916 


(10) 


20 


170.363 


1.033 


(20) 


170.363 


1.033 


(20) 


170.026 


1.033 


(20) 


30 


380.045 


1.104 


(30) 


380.045 


1.104 


(30) 


379.611 


1.104 


(30) 


40 


666.975 


1.156 


(40) 


666.975 


1.156 


(40) 


666.467 


1.156 


(40) 


50 


1028.596 


1.197 


(50) 


1028.596 


1.197 


(50) 


1027.998 


1.197 


(50) 


60 


1462.912 


1.231 


(60) 


1462.912 


1.231 


(60) 


1462.210 


1.231 


(60) 


70 


1968.278 


1.261 


(70) 


1968.278 


1.261 


(70) 


1967.484 


1.261 


(70) 


80 


2543.311 


1.287 


(80) 


2543.311 


1.287 


(80) 


2542.467 


1.287 


(80) 


90 


3186.982 


1.310 


(90) 


3186.983 


1.310 


(90) 


3185.995 


1.310 


(90) 


100 


3898.102 


1.331 


(100) 


3898.103 


1.331 


(100) 


3897.051 


1.331 


(100) 


110 


4675.652 


1.355 


(2,108) 


4675.858 


1.355 


(2,108) 


4674.619 


1.355 


(2,108) 


120 


5518.467 


1.375 


(4,116) 


5518.839 


1.376 


(4,116) 


5517.415 


1.374 


(3,117) 


130 


6425.708 


1.395 


(6,124) 


6426.152 


1.396 


(6,124) 


6424.708 


1.396 


(6,124) 


140 


7396.731 


1.413 


(8,132) 


7397.271 


1.414 


(8,132) 


7395.796 


1.414 


(8,132) 


150 


8430.934 


1.433 


(12,138) 


8431.678 


1.433 


(12,138) 


8430.054 


1.432 


(11,139) 


160 


9527.713 


1.448 


(14,146) 


9528.564 


1.449 


(14,146) 


9526.851 


1.447 


(13,147) 


170 


10686.607 


1.464 


(17,153) 


10687.544 


1.464 


(17,153) 


10685.664 


1.463 


(16,154) 


180 


11906.937 


1.479 


(20,160) 


11907.934 


1.477 


(18,162) 


11905.975 


1.478 


(19,161) 


190 


13188.284 


1.492 


(22,168) 


13189.315 


1.492 


(22,168) 


13187.303 


1.492 


(22,168) 


200 


14530.226 


1.506 


(26,174) 


14531.342 


1.504 


(24,176) 


14529.198 


1.506 


(25,175) 



TABLE I: Lowest energy En, and outer shell radius Rn found for ion clusters in the octupole trap, as obtained from Monte 
Carlo global minimization (left columns), from the discrete (Thomson) or continuous (asymptotic) shell models (central and 
right columns, respectively). The ion arrangements into shells predicted by the three methods are also indicated. 



as a function of the number of ions in Fig. [TJ The av- 
erage radii obtained from globally optimized structures, 
superimposed for selected sizes N — I0 k , k = 1-4, cor- 
rectly match Eq. (J6]) until about 1000 ions are reached, 
larger clouds showing clear deviations in the inner shells. 
However, the cloud radius (outermost shell) and the min- 
imum energy are both accurately reproduced. One pre- 
diction of the continuous shell model is the size at which 
new shells appear. In contrast with the quadrupole case, 
where shells are essentially added over an existing core 
[Til [l9| , shells for octupole clusters grow both on the out- 
side and on the inside. The shell capacitance may then 
be defined as the maximum number of ions that a shell 
can sustain, and above which a new shell appears. The 
continuous shell model predicts the onset of new shells 
at N = 109, N = 442, N = 1129, N = 2264, N = 3992, 
N = 6466, N = 9709, N = 13967, and N = 19249. The 
model also predicts that the 10 5 -ion cluster should have 
17 shells but, as will be seen below, structural minimiza- 
tion yields a more contrasted picture. These transitions 
occur at much larger sizes than in the quadrupole case 
[l8l Qjij], where for instance the 13- and 58-ion clusters 
adopt the (1,12) and (1,12,45) shell structures, respec- 
tively [13|. 

Up to 108 ions, the most stable configurations are the 
same as the Thomson minima (29l . [30| , within some mi- 
nor radial relaxations. The Thomson shell model also 
performs rather well for predicting structures of clusters 
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FIG. 1: (Color online) Shell radii versus cluster size, as pre- 
dicted from the correlated shell model (solid lines; different 
curves correspond to successive shells) and from global op- 
timization (open circles). The largest shell radius follows a 
iV 1 / 7 scaling law. Inset: lowest energy versus cluster size, as 
obtained from the correlated shell model (dashed line) and 
from global optimization (open circles). 



containing two shells. This agreement suggests to use 
the known Thomson minima for optimizing the atomic 
structure of multishell systems as well, as long as all shells 
remain very thin. Such an idea is also supported by a re- 
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cent investigation by Cioslowski on a modified but sim- 
ilar Thomson model [33j. We have guided Monte Carlo 
global minimization by combining the results of the shell 
model with the available solutions to the Thomson prob- 
lem [301 ] . In this approach, the continuous shell model 
is used to predict the optimal number of layers and the 
individual numbers of ions per layer (the discrete shell 
model could also be used). The coordinates of the ions 
in each layer are then scaled from the corresponding so- 
lution to the Thomson problem to have a radius given by 
the predictions of the shell model, and the only remaining 
degrees of freedom are two Euler angles {9i, fa} for each 
layer. The new Monte Carlo optimization consists thus 
in first locating the low-energy regions in these angles 
space, and to locally minimize the resulting structures 
by relaxing all ionic positions. This method was found 
to perform very satisfactorily with respect to brute force 
basin-hopping minimization for several sizes in the range 
N = 100-1000, leading to low-energy structures often 
identical, or higher but by only a few 10 -3 percents. 

Provided that correlations are included in the evalua- 
tion of £c? tra (jV), the above results show that the com- 
bined shcll+Montc Carlo optimization should be espe- 
cially useful for large systems, for which successful basin- 
hopping runs would be prohibitive, and this even allowed 
us to explore clusters containing up to N = 10 5 ions. At 
large sizes, Fig. [1] shows that the cloud radius R and its 
minimum energy $ scale as R oc jV 1 / 7 and $ cx TV 13 / 7 , 
respectively. Both scaling laws readily follow from a sim- 
ple cold fluid approximation, which is expected to be 
valid for large sizes [H, [35|. The maximum radius R 
of the cloud naturally relates to the total number N of 
ions through a n(r) oc r A radial density. For the energy 
$, the Coulomb and trapping components are related 
to each other through the virial theorem as <&c = 6$t 
at any local minimum configuration where the gradient 
vanishes. Hence $ is proportional to the trapping energy 
only, which is simply integrated over the radial density 
n(r) as $ = 7$ T oc R 13 oc N 13 / 7 . 

The validity of the cold fluid theory is better mani- 
fested on the radial effective density profile. From the 
stable ionic configurations obtained for the sizes N = 
10 fc , k = 3-5, we have calculated the minimum pair 
distance y™ m between a given ion i and all other ions 
j, which is related to the local density n(r) through 
r™j ln ~ n^ 1 / 3 . Fig. [2] shows the correlation between 
this quantity and the radial distance Ti of ion i for the 
three sizes. The r -4 / 3 dependence highlighted in Fig. [2] 
thus confirms that the cold fluid theory holds increasingly 
well for these ion clouds, deviations being more notice- 
able in the less populated (but more fluctuating) inner 
regions. The most striking feature of Fig. [5] is the broad- 
ening of the shells when the number of ions increases 
from 10 3 to 10 4 , especially noticeable for the inner shells. 
Above a few thousands of ions (an estimate based on 
mixed shell/Thomson optimization yields N* ~ 4000), 
this broadening is sufficient for some shells to overlap 
into a more continuous radial distribution of ions below 
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FIG. 2: (Color online) Minimum pair distance versus radial 
distance in stable structures containing 10 3 -10 5 ions. Inset: 
fraction of ions inside a sphere of given radius. The most 
stable configuration of the 1000-ion cluster is also depicted, 
with the front quarter removed. 

the outermost layer, which is the only one to remain thin, 
even for the largest size considered, N — 10 5 . This mixed 
continuous/discrete behavior is better seen on the accu- 
mulated fraction x( r ) of ions inside a sphere of radius 
r. The variations of x with increasing r, depicted as an 
inset in Fig. [2J exhibit spectacular changes as the size 
reaches 10 4 . The sharp steps found in small clouds or 
at large radii are characteristic of new shells, but are 
progressively softened in the inner regions of the 10 4 -ion 
system and even replaced by a nearly continuous profile 
in the largest cluster. In an octupole trap, crystallization 
should then be understood as the formation of an outer- 
most thin layer with a softer, decreasingly dense but thick 
inner layer. Surface effects, which play a major role in 
reducing the melting point of ion clouds in harmonic po- 
tentials [22| , could thus have a very different influence in 
the case of the octupole trap. 

III. FINITE TEMPERATURE PROPERTIES 

Whereas the previous section considered the stable 
structures and static issues, we now discuss the thermo- 
dynamical and dynamical behavior of selected ion clus- 
ters in octupole traps. Classical MD simulations have 
been carried out at various energies to compute sev- 
eral thermodynamical and structural observables. Or- 
der/disorder transitions have been monitored using the 
root mean square bond length fluctuation or Lindemann 
index 5, as well as the particle-resolved diffusion constant 
Di from the integrated velocity autocorrelation function. 
The rms bond length fluctuation index has been com- 
puted for collective parts of the clouds (intra-shell in- 
dices) or from pairs of ions belonging to different shells 
in the stable structure (inter-shell index) . The variations 
of the shell-resolved Lindemann indices with the kinetic 
temperature are shown in Fig. [3] for the 512-ion system 
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made of 3 shells. This size is small enough for the inner 
shells to be well defined and not overlapping with each 
other. All Lindemann indices exhibit smooth variations 
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FIG. 3: (Color online) Relative bond length fluctuation in- 
dices obtained for the 512-ion cluster from molecular dynam- 
ics simulations. The dotted vertical line marks the onset of 
the melting temperature for this cluster (T m ~ 1.1 x 10~ 2 ). 
Inset: dependence of the melting temperature with inverse 
cluster radius, and extrapolation to the infinitely large limit 
based on results obtained for N > 512. 



at low and high temperatures, and one sharp increase 
above S > 0.15 in a narrow temperature range, allowing 
an accurate estimation of the corresponding melting tem- 
perature. From Fig. [31 the global melting of the 512-ion 
system can therefore be located near T m ~ 1.1 x 10~ 2 . 
However, the additional Lindemann indices clearly reveal 
that the interior shell exhibits some preliminary soften- 
ing already above T ~ 10~ 3 before fully melting near 
Tm — 3 x 10~ 3 . Melting of this inner shell impacts 
the distance fluctuations within the other shells, as seen 
from the slight increase in the corresponding Lindemann 
indices. However, these shells clearly undergo their own 
distinct melting transitions at Tm — 8.5 x 10~ 3 (second 
shell) and t! 3) ~ 1.4 x 1CT 2 (outermost shell). The in- 
tershell Lindemann index follows similar variations as the 
intrashell index of the intermediate, second shell. This 
agreement is not fortuitous: melting of the second shell 
occurs while the interior shell is already disordered, and 
the present results show that these two shells tend to 
mix, to some extent, while the outermost layer remains 
rigid. 

The melting mechanisms can be analysed in more de- 
tails by looking at the ions motion at various tempera- 
tures. The correlations between the average radial dis- 
tance (rj) and the diffusion constant Di are shown in 
Fig. 0] for the 2048-ion system at four characteristic tem- 
peratures. For this system, the melting point obtained 
from the variations of the global Lindemann index lies 
near T m ~ 0.023. At T — 0.011, the four shells are 
clearly seen as narrow vertical spots with low values for 
Di at radii close to 1.32, 1.68, 1.97, and 2.20. The inter- 



shell motion and occasional hops of ions between neigh- 
boring shells are seen at T — 0.018, they are associated 
with a much higher diffusion constant. Note that Di 
exhibit a steady decrease with the radial distance, in 
agreement with the previously found softening of inter- 
nal layers. At T — 0.028 the global diffusivity is also 
higher, and the shell structure seems essentially lost. A 
similar trend is found at the highest temperature consid- 
ered here, T = 0.048, with a single broad spot centered 
around a radius of 2.02. The peculiar core melting ef- 
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FIG. 4: (Color online) Ion diffusion constant as a function of 
its time average radius, in the 2048-ion cluster at four tem- 
peratures. 

feet is a consequence of the lower density of the inner 
shells, and is similar in this respect to the general sur- 
face melting process in solid state and nanoscale mate- 
rials [3o| . At first sight, this phenomenon may preclude 
from unambiguously defining the melting temperature of 
the entire system since, strictly speaking, the system is 
not yet fully melted when the global Lindemann index 
barely exceeds 0.15. However, because the relative num- 
ber of ions in the external regions grows with size as 
N 1 / 1 according to the cold fluid model, the definition of 
the global index 5 should reflect more and more closely 
the value of the outermost layer. The melting point T m 
was thus defined for all sizes as the temperature at which 
<5 exceeds 0.15. For the present clouds in octupole traps, 
the variations of T m have been represented in the inset of 
Fig. [3] as a function of inverse cluster radius. Not surpris- 
ingly, fluctuations are seen at small sizes, and the 32-ion 
system (1/i? ~ 0.61) seems extra resistant to melting 
relative to its neighboring sizes. Since the 31- and 33- 
ion systems have a much lower melting point, the special 
stability of the 32-ion cluster indicates a magic charac- 
ter, further confirmed by its high icosahedral symmetry. 
As N reaches 256 (1/i? — 0.45), T m increases linearly 
as the inverse radius decreases. This linear depression, 
though not anticipated for these highly heterogeneous, 
core-melted systems, allows some straightforward extrap- 
olating to the infinite size limit 1/i? — > 0, leading to 
Too — 0.076. Even though it is hard to figure out what 
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a cloud confined in an octupole trap would physically 
represent at the bulk limit, the present simulation re- 
sults should be very useful for estimating melting tem- 
peratures in system sizes of the order of 10 5 or more, as 
experimentally studied in other groups [ill ]. 

IV. CONCLUDING REMARKS 

Despite their obvious differences, ion clouds in oc- 
tupole and quadrupole confinements share several re- 
markable features. Firstly, the onset of the transition 
from a multishcll structure to the bcc Wigner crystal in 
quadrupole trap was estimated to be around 10 4 ions p^ |. 
In octupole traps the shell structure becomes blurred 
above a few thousands of ions. Secondly, for both cases, 
the depression in the melting point with respect to the 
bulk limit follows a linear scaling with inverse cluster 
radius [22]. Thirdly, in small systems, a clear differen- 
tial melting between the inner and outer shells seems to 
take place in both confinements. These similarities may 
well extend past the specific octupole trap, hold also in 
higher-order traps, and even be universal. As the expo- 
nent p = 6 in Eq. takes higher integer values, the sta- 
ble shell structure should be qualitatively preserved, only 
with fewer shells. Therefore the transition to a softened 
core should be delayed as p increases, but one cannot 
exclude that the same melting mechanisms will remain, 
and in particular that the melting temperature will dis- 
play linear variations with inverse radius. 

Because ion clouds in octupole traps remain poorly 
studied, the present work could be extended along many 
lines. Keeping the isotropic case as a model of a more re- 



alistic three-dimensional trap, it would be interesting to 
examine more specifically dynamical properties involving 
the vibrations and normal modes [33, [38| or the transi- 
tions between regular and chaotic regimes [H, H(| ■ Ex- 
tension of analytical models [3l|, [H| developed for the 
quadrupole case to higher-order traps would also be use- 
ful. 

Finally, at the price of introducing additional param- 
eters in the model, one natural step beyond the present 
work would be to look at clusters in linear octupole traps, 
such as those discussed by Okada and coworkers [24| ■ In 
these systems, confinement along the symmetry axis is 
harmonic in nature, whereas confinement perpendicular 
to the axis is octupolar. The stable structures are hollow 
tubes of ions [24| , and it is unclear how the combined con- 
tributions of the harmonic and octupolar confinements 
will determine the cluster properties. Even aware of 
these complications, the multishell structure and core- 
melted phase predicted here, as well as the scaling laws 
connecting size, radius, and melting temperature, should 
all become amenable to experimental comparison in the 
near future. May such measurements also help us under- 
standing the collective properties of these exotic states of 
Coulombic matter. 
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